#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>

static void _rotating_calipers(const std::vector<cv::Point2f>& hull, cv::Point2f out[3]) {
    int n = hull.size();
    CV_Assert(n > 2);
    
    int left = 0, bottom = 0, right = 0, top = 0;
    float left_x, right_x, top_y, bottom_y;
    cv::Point2f pt0 = hull[0];
    left_x = right_x = pt0.x;
    top_y = bottom_y = pt0.y;
    
    std::vector<cv::Point2f> vect(n);
    std::vector<float> invVect(n);
    for (int i = 0; i < n; i++) {
        if (pt0.x < left_x) {
            left_x = pt0.x;
            left = i;
        }
        if (pt0.x > right_x) {
            right_x = pt0.x;
            right = i;
        }
        if (pt0.y < bottom_y) {
            bottom_y = pt0.y;
            bottom = i;
        }
        if (pt0.y > top_y) {
            top_y = pt0.y;
            top = i;
        }
        
        cv::Point2f pt = hull[(i+1 >= n ? 0 : i+1)];
        double dx = pt.x - pt0.x;
        double dy = pt.y - pt0.y;
        vect[i].x = (float)dx;
        vect[i].y = (float)dy;
        invVect[i] = (float)(1./std::sqrt(dx*dx+dy*dy));
        
        pt0 = pt;
    }
    
    float direction = 0;
    double ax = vect[n-1].x;
    double ay = vect[n-1].y;
    for (int i = 0; i < n; i++) {
        double bx = vect[i].x;
        double by = vect[i].y;
        double convexity = ax*by-ay*bx;
        if (convexity != 0) {
            direction = (convexity > 0) ? 1.f : -1.f;
            break;
        }
        ax = bx;
        ay = by;
    }
    
    int seq[4];
    seq[0] = bottom;
    seq[1] = right;
    seq[2] = top;
    seq[3] = left;
    
    float a = direction;
    float b = 0;
    
    float minArea = FLT_MAX;
    float minAreaA, minAreaB, minAreaWidth, minAreaHeight;
    int minAreaLeft, minAreaBottom;
    
    for (int k = 0; k < n; k++) {
        float dotProduct[4];
        dotProduct[0] = a*vect[seq[0]].x + b*vect[seq[0]].y;
        dotProduct[1] = -b*vect[seq[1]].x + a*vect[seq[1]].y;
        dotProduct[2] = -a*vect[seq[2]].x - b*vect[seq[2]].y;
        dotProduct[3] = b*vect[seq[3]].x - a*vect[seq[3]].y;
        
        float maxcos = dotProduct[0] * invVect[seq[0]];
        int mainElement = 0;
        for (int i = 1; i < 4; i++) {
            float temp = dotProduct[i] * invVect[seq[i]];
            if (temp > maxcos) {
                maxcos = temp;
                mainElement = i;
            }
        }
        
        int idx = seq[mainElement];
        float xx = vect[idx].x * invVect[idx];
        float yy = vect[idx].y * invVect[idx];
        switch (mainElement) {
            case 0:
                a = xx;
                b = yy;
                break;
            case 1:
                a = yy;
                b = -xx;
                break;
            case 2:
                a = -xx;
                b = -yy;
                break;
            case 3:
                a = -yy;
                b = xx;
                break;
        }
        
        seq[mainElement] += 1;
        if (seq[mainElement] == n) {
            seq[mainElement] = 0;
        }
        
        float dx = hull[seq[1]].x - hull[seq[3]].x;
        float dy = hull[seq[1]].y - hull[seq[3]].y;
        float width = dx*a + dy*b;
        dx = hull[seq[2]].x - hull[seq[0]].x;
        dy = hull[seq[2]].y - hull[seq[0]].y;
        float height = -dx*b + dy*a;
        float area = width * height;
        if (area <= minArea) {
            minArea = area;
            minAreaA = a;
            minAreaB = b;
            minAreaWidth = width;
            minAreaHeight = height;
            minAreaLeft = seq[3];
            minAreaBottom = seq[0];
        }
    }
    
    float a1 = minAreaA;
    float b1 = minAreaB;
    float a2 = -b1;
    float b2 = a1;
    float c1 = a1*hull[minAreaLeft].x + hull[minAreaLeft].y * b1;
    float c2 = a2 * hull[minAreaBottom].x + hull[minAreaBottom].y*b2;
    float idet = 1.f / (a1*b2-a2*b1);
    float px = (c1*b2-c2*b1)*idet;
    float py = (a1*c2-a2*c1)*idet;
    out[0].x = px;
    out[0].y = py;
    out[1].x = a1*minAreaWidth;
    out[1].y = b1*minAreaWidth;
    out[2].x = a2*minAreaHeight;
    out[2].y = b2*minAreaHeight;
}

static cv::RotatedRect _min_area_rectangle(const std::vector<cv::Point>& pts) {
    std::vector<cv::Point> temp;
    cv::convexHull(pts, temp, true);
    std::vector<cv::Point2f> hull;
    for (int i = 0; i < temp.size(); i++) {
        hull.push_back(cv::Point2f(temp[i].x, temp[i].y));
    }
    int n = hull.size();
    cv::RotatedRect rotatedRect;
    if (n > 2) {
        cv::Point2f out[3];
        _rotating_calipers(hull, out);
        
        rotatedRect.center.x = out[0].x + (out[1].x + out[2].x) * 0.5f;
        rotatedRect.center.y = out[0].y + (out[1].y + out[2].y) * 0.5f;
        rotatedRect.size.width = (float)std::sqrt((double)out[1].x * out[1].x + (double)out[1].y*out[1].y);
        rotatedRect.size.height = (float)std::sqrt((double)out[2].x*out[2].x + (double)out[2].y*out[2].y);
        rotatedRect.angle = (float)std::atan2((double)out[1].y, (double)out[1].x);
        
    } else if (n == 2) {
        rotatedRect.center.x = (hull[0].x + hull[1].x) * 0.5f;
        rotatedRect.center.y = (hull[0].y + hull[1].y) * 0.5f;
        double dx = hull[1].x - hull[0].x;
        double dy = hull[1].y - hull[0].y;
        rotatedRect.size.width = (float)std::sqrt(dx*dx+dy*dy);
        rotatedRect.size.height = 0;
        rotatedRect.angle = (float)std::atan2(dy, dx);
    } else if (n == 1) {
        rotatedRect.center = hull[0];
    }
    rotatedRect.angle = (float)rotatedRect.angle* 180 / CV_PI;
    
    return rotatedRect;
}

int main(int argc, char **argv) {
    cv::Mat src = cv::imread("/home/xlll/Downloads/opencv/samples/data/lena.jpg", cv::IMREAD_COLOR);
    if (src.empty()) {
        std::cout << "failed to read image!" << std::endl;
        return EXIT_FAILURE;
    }
    
    cv::Mat binary;
    //cv::Canny(src, binary, 250, 1000);
    cv::Canny(src, binary, 250, 600);
    
    std::vector<std::vector<cv::Point>> contours;
    cv::findContours(binary, contours, cv::RETR_LIST, cv::CHAIN_APPROX_NONE);
    std::vector<cv::RotatedRect> rotatedRects, rotatedRects1;
    std::vector<cv::Rect> rects;
    for (int i = 0; i < contours.size(); i++) {
        rotatedRects.push_back(cv::minAreaRect(contours[i]));
        rotatedRects1.push_back(_min_area_rectangle(contours[i]));
        rects.push_back(cv::boundingRect(contours[i]));
        for (int j = 0; j < contours[i].size(); j++) {
            cv::circle(src, contours[i][j], 1, cv::Scalar(0, 255, 255), cv::FILLED);
        }
    }
    
    bool allSame = true;
    for (int i = 0; i < rotatedRects.size(); i++) {
        cv::Point2f pts[4];
        rotatedRects[i].points(pts);
        std::cout << "No. " << i << ": angle = " << rotatedRects[i].angle << ", center = " << rotatedRects[i].center << ", size = " << rotatedRects[i].size << std::endl;
        for (int j = 0; j < 4; j++) {
            cv::line(src, pts[j], pts[(j+1)%4], cv::Scalar(255, 0, 0), 2, cv::LINE_AA);
        }
        cv::rectangle(src, rects[i], cv::Scalar(0, 255, 0), 2, cv::LINE_AA);
        
        if (allSame) {
            if (std::fabs(rotatedRects[i].angle - rotatedRects1[i].angle) > FLT_EPSILON ||
                std::fabs(rotatedRects[i].center.x - rotatedRects1[i].center.x) > FLT_EPSILON ||
                std::fabs(rotatedRects[i].center.y - rotatedRects1[i].center.y) > FLT_EPSILON ||
                std::fabs(rotatedRects[i].size.width - rotatedRects1[i].size.width) > FLT_EPSILON ||
                std::fabs(rotatedRects[i].size.height - rotatedRects1[i].size.height) > FLT_EPSILON) {
                allSame = false;
            }
        }
    }
    std::cout << "all results are " << (allSame ? "the same!" : "not the same!") << std::endl;
    
    cv::imshow("src", src);
    cv::imshow("binary", binary);
    cv::waitKey(0);
    
    return 0;
}
